Comment modéliser une relation linéaire Y = a.X+b ?
Claude Grasland & Jean-Paul Nguesso
2025-06-01
1. PREPARATION DES DONNEES
1.1. Chargement du tableau principal
On charge le fichier des pays en 2018
don <-read.table(file ="data/africa_pays_2018/data/africa_pays_2018.csv", # nom du fichier et chemin d'accèssep =";", # séparateur (ici, des points-virgule)dec=",", # Type de décimaleheader =TRUE, # ligne d'en-tête avec le nom des variablesencoding="UTF-8") # encodage adapté au français
1.2 Choix des deux variables à analyser
En dehors de iso3 et nom, on ne garde que deux variables que l’on renomme X et Y avec colnames() et que l’on convertit en type numérique général. Il suffira par la suite de modifier le choix des variables X et Y pour faire d’autres analyses.
CODE NOM X Y
1 AGO Angola 65.9 14.3
2 BDI Burundi 13.2 2.7
3 BEN Bénin 47.6 20.0
4 BFA Burkina Faso 29.7 16.0
5 BWA Botswana 69.8 47.0
6 CAF Rep. Centrafricaine 41.6 4.3
#par(mfrow=c(1,1))#monplot(sel$X,sel$Y,sel$CODE)
1.3 On est malin …
Mais comme on ne sait plus ce que sont X et Y, on le précise avec des chaînes de caractères qu’on pourra utiliser dans les graphiques. Et on peut préparer une version multilangue …
# Pour la version françaisetitre <-"Les pays d'Afrique en 2018"nomX <-"Taux d'urbanisation (%)"nomY <-"Taux d'accès à internet (%)"auteur <-"Ecole d'été AfromapR, Bouaké 2025"
2. NORMALITE, VISUALISATION, CORRELATION
2.1 Vérification de la normalité de X et Y
La régression linéaire met en relation deux variables quantitatives X et Y dont on suppose que la distribution est normale (gaussienne) , c’est-à-dire unimodale et symérique.
Une analyse rapide avec boxplot() permet de repérer s’il y a des valeurs exceptionnelles.
Pour X comme Y, on ne voit pas apparaître de valeurs exceptionnelles c’est-à-dire de points situés au delà de l’intervalle défini par la boîte à moustache.
2.1 Vérification de la normalité de X et Y
Les fonctions qqnorm() et qqline() sont plus précises …
Y ne suit pas bien la droite donc elle n’est pas tout à fait gaussienne.
2.1 Vérification de la normalité de X et Y
Mais la solution la plus précise est le test de Shapiro qui pose l’hypothèse H0 : la distribution est normale.
shapiro.test(sel$X)
Shapiro-Wilk normality test
data: sel$X
W = 0.97546, p-value = 0.3931
shapiro.test(sel$Y)
Shapiro-Wilk normality test
data: sel$Y
W = 0.89916, p-value = 0.0005157
Si le test est supérieur à 0.05, on peut considérer que la distribution est gaussienne ce qui est le cas de X. Dans le cas contraire (p <0.05) la distribution n’est pas gaussienne ce qui est le cas de Y.
Mais en pratique on accepte en général de faire des analyses de régression avec des distributions non-gaussienne, surtout s’il n’y a pas de valeurs exceptionnelles.
2.2 Visualisation de la forme de la relation
On peut faire un simple plot(X,Y). Mais on peut aussi créer pour cela une fonction personalisée adapté à ses préférences
monplot <-function (varX , varY, varN ){ Ymin <-min(varY) Ymax <-max(varY) Ysup <- Ymax+(Ymax-Ymin)*0.1plot(varX,varY,main = titre, # titrecex.main =1, # police du titrecex =0.6, # taille des symbolespch =19, # cercles pleinsylim =c(Ymin, Ysup),col ="red") # couleur des symbolestext(varX,varY,varN,cex=0.5,pos=3) # nom des élémentabline(v=mean(varX),lty=2,lwd=1,col="blue") # moyenne Xabline(h=mean(varY),lty=2,lwd=1,col="blue") # moyenne Ygrid() }
2.2 Visualisation de la forme de la relation
Je peux désormais utiliser ma fonction monplot() !
Je peux décider de ne pas afficher le label des points.
monplot(varX = sel$X,varY = sel$Y, varN =NULL)
2.3 Analyse de la corrélation
Je commence par calculer le coefficient de corrélation linéaire (r) et le pouvoir explicatif de X par rapport à Y (r2)
cor(sel$X,sel$Y) # coefficient de corrélation (r)
[1] 0.530947
100*cor(sel$X,sel$Y)**2# pouvoir explicatif (r2)
[1] 28.19048
2.3 Analyse de la corrélation
Puis, je teste la significativité de la corrélation linéaire …
cor.test(sel$X,sel$Y) # test de significativité (p-value)
Pearson's product-moment correlation
data: sel$X and sel$Y
t = 4.2955, df = 47, p-value = 8.68e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.2935825 0.7066417
sample estimates:
cor
0.530947
2.3 Analyse de la corrélation
… et je la compare à celle du coefficient de corrélation de rang de Spearman
cor.test(sel$X,sel$Y, method="spearman") # test de significativité (p-value)
Spearman's rank correlation rho
data: sel$X and sel$Y
S = 10474, p-value = 0.0007479
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.4656206
Les deux coefficients sont assez proches (+0.53 pour Pearson et +0.47 pour Spearman) ce qui est rassurant. En général, si les deux coefficients diffèrent il y a un problème de valeur exceptionnelle ou de relation non linéaire.
2.3 Analyse de la corrélation
On peut conclure des analyses précédentes que :
il existe une relation significative (p-value < 0.01)
cette relation est positive (r = +0.53 )
cette relation a un pouvoir explicatif assez faible (r2 = 28%)
la relation est monotone et globalement linéaire car le coefficient de Spearman est à peu près égale au coefficient de Pearson. Il n’y a donc pas besoin de transformer les variables ou retirer des valeurs exceptionnelles.
3. LA REGRESSION LINEAIRE
3.1 Hypothèses statistiques
Conditions a priori
X et Y sont deux variables normales (gaussienne)
il existe une corrélation significative entre X et Y (p< 0.05)
X explique une part suffisamment forte de Y (r2 > 20% )
Le nuage de point affiche une forme linéaire
les points sont répartis de façon régulière le long du nuage de points
Il n’y a pas de valeurs exceptionnelles susceptibles de perturber le calcul.
3.1 Hypothèses statistiques
Méthode des moindres carrés ordinaire (MCO)
La droite \(y_i = a.x_i + b + \epsilon_i\) qui minimise la somme des carrés des écarts entre les valeurs observées \(y_i\) et les valeurs estimées \(\hat{y_i}\) a pour équation :
La somme des carré des écarts totale (\(SCE_{tot}\)) correspond à la variance de la variable à expliquer : \(SCE_{tot} = \sum_{i=1}^k (y_{i}-\bar{y})^2\)
La somme des carré des écarts résiduels (\(SCE_{err}\)) correspond à la somme des carrés des différences entre valeurs observées et estimées \(SCE_{error} = \sum_{i=1}^k (y_{i}-\hat{y})^2\)
Le pouvoir explicatif d’un modèle de régression correspond à la part de la variance de Y expliquée par X.
Un modèle de régression n’est valide que si les résidus de ce modèle \(\epsilon_i = (y_i - \hat{y}_i)\) remplissent les conditions suivantes :
Normalité de la distribution des résidus
Absence d’autocorrélation des résidus
Homogénéité de la variance des résidus
Absence de valeur à fort effet de levier
Si ces quatre conditions ne sont pas remplies, les estimations de Y en fonction de X seront entâchées d’erreur et leur intervalle de confiance ne sera pas valable.
3.2 La fonction lm()
La fonction lm() ou lm est l’abbréviation de linear model permet d’effectuer la plupart des modèles de régression linéaire basés sur la méthode des moindres carrés ordinaire. Sa syntaxe est a priori très simple et renvoie les coefficients b et a du modèle de régression.
Cette droite a pour équation Y = 0.5X + 2.8 ce qui veut dire que chaque fois que l’urbanisation augemente de 1 pt l’usage d’internet augmente de 0.5 pts. Un pays qui aurait une urbanisation de 50% devrait donc avoir un taux d’usage d’internet égal à 0.5 x 50 + 2.8 = 27.8%
3.5 Analyse des résidus
On peut extraire de l’objet créé par lm() les valeurs estimées de Y et les résidus c’est-à-dire les erreurs d’estimation.
sel$Y_estim<-monmodel$fitted.values # Valeurs estiméessel$Y_resid<-monmodel$residuals # Valeurs résiduellessel <- sel[order(sel$Y_resid),] # Tri du tableau en fonction des résidussel
Les résidus négatifs correspondent aux pays qui ont moins accès à internet que ce que leur urbanisation laisserait prévoir :
head(sel,5)
CODE NOM X Y Y_estim Y_resid
10 COG Congo 67.2 8.7 37.03663 -28.33663
40 SOM Somalie 45.3 2.0 25.33241 -23.33241
24 LBY Libye 80.3 21.8 44.03778 -22.23778
1 AGO Angola 65.9 14.3 36.34186 -22.04186
14 ERI Erythrée 40.4 1.3 22.71366 -21.41366
Ainsi le Congo n’a que 8.7% d’accès à internet (Y) alors que son taux d’urbanisation est de 67.2% (X) ce qui laisserait prévoir un taux d’accès à internet de 37.0% (Y_estim). Il a donc un résidus de -28 pts (Y_resid) qui signifie que son taux d’accès à internet est beaucoup plus faible que ce que laisserait prévoir son urbanisation.
3.5 Analyse des résidus
Les résidus positifs correspondent aux pays qui ont plus accès à internet que ce que leur urbanisation laisserait prévoir :
tail(sel,5)
CODE NOM X Y Y_estim Y_resid
32 NAM Namibie 50.5 51.0 28.11149 22.88851
13 EGY Egypte 42.7 46.9 23.94287 22.95713
44 TUN Tunisie 69.1 64.2 38.05206 26.14794
26 MAR Maroc 62.8 64.8 34.68509 30.11491
41 SWZ Swaziland 23.9 47.0 13.89542 33.10458
-Ainsi le Swaziland a 47 % d’accès à internet (Y) alors que son taux d’urbanisation est de 23.9% (X) ce qui laisserait prévoir un taux d’accès à internet de 13.9% (Y_estim). Il a donc un résidus de +33 pts (Y_resid) qui signifie que son taux d’accès à internet est beaucoup plus fort que ce que laisserait prévoir son urbanisation.
4. COMPLEMENTS STATISTIQUES (facultatifs)
4.1 Hypothèses statistiques
Vérifications a posteriori
Un modèle de régression n’est valide que si les résidus de ce modèle \(\epsilon_i = (y_i - \hat{y}_i)\) remplissent les conditions suivantes :
Normalité de la distribution des résidus
Absence d’autocorrélation des résidus
Homogénéité de la variance des résidus
Absence de valeur à fort effet de levier
Si ces quatre conditions ne sont pas remplies, les estimations de Y en fonction de X seront entâchées d’erreur et leur intervalle de confiance ne sera pas valable.
On charge le package car (companion to applied regession) pour réaliser certains des diagnostics
library(car)
4.1 Diagnostic 1 : Indépendance des résidus ?
L’objectif est de savoir si les résidus se répartissent régulièrement de part et d’autre de la droite de régression tout au long de celle-ci. Si c’est bien le cas le graphique residuals Vs Fitted généré par plot(monmodel,1) devrait donner une droite horizontale :
plot(monmodel,1,labels.id = sel$CODE)
4.1 Diagnostic 1 : Indépendance des résidus ?
4.1 Diagnostic 1 : Indépendance des résidus ?
On peut tester statistiquement l’indépendance des résidus à l’aide du test de Durbin-Watson qui mesure si deux valeurs successives ont des résidus proches. L’indépendance des résidus est rejetée si p-value < 0.05
durbinWatsonTest(monmodel)
lag Autocorrelation D-W Statistic p-value
1 0.00312332 1.94811 0.856
Alternative hypothesis: rho != 0
Ici on trouve p-value > 0.05 donc les résidus sont indépendants.
Diagnostic 2 : Normalité des résidus ?
L’objectif est de savoir si les résidus ont une distribution normale Si c’est bien le cas le graphique généré par plot(monmodel,2) devrait donner une droite oblique :
plot(monmodel,2,labels.id = sel$CODE)
Diagnostic 2 : Normalité des résidus ?
Diagnostic 2 : Normalité des résidus ?
On peut tester statistiquement la normalité des résidus à l’aide du test de Shapiro-Wilk. Les résidus sont normaux si p-value > 0.05
shapiro.test(monmodel$residuals)
Shapiro-Wilk normality test
data: monmodel$residuals
W = 0.96515, p-value = 0.1543
Ici on trouve une p-value supérieure à 0.05 donc la distribution des résidus est approximativement gaussienne.
Diagnostic 3 : Homogénéité des résidus ?
L’objectif est de savoir si la variance des résidus est constante, c’est-à-dire si il s’écarte environ de la même distance tout au long de la droite . Si c’est bien le cas le graphique généré par plot(monmodel,3) devrait donner une droite horizontale
plot(monmodel,3,labels.id = sel$CODE)
Diagnostic 3 : Homogénéité des résidus ?
Diagnostic 3 : Homogénéité des résidus ?
On peut tester statistiquement l’homogénéité des résidus à l’aide du test de Breush-Pagan. L’hypothèse d’homogénéité est rejetée si la p-value est inférieure à 0.05.
ncvTest(monmodel)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 3.029338, Df = 1, p = 0.081771
Ici, la p-value est supérieure à 0.05 donc les résidus sont globalement homogènes.
Diagnostic 4 : Absence de valeur exceptionnelles ?
L’objectif est de savoir s’il existe des valeurs qui exercent une influence exceptionnelle sur les résultats de la régression. On peut reprérer ces valeurs de plusieurs manières, notamment à l’aide de la distance de Cook générée par plot(monmodel,4).
plot(monmodel,4,labels.id = sel$CODE)
Diagnostic 4 : Absence de valeur exceptionnelles ?
On repère les pays les plus susceptibles de perturber le modèle (Sierra Léone, Sud Soudan et Côte d’Ivoire)
Diagnostic 4 : Absence de valeur exceptionnelles ?
Le test statistique de Bonferroni permet de déterminer s’il existe des valeurs exceptionnelles avec une p-value < 0.05.
outlierTest(monmodel, labels = sel$CODE)
No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferroni p
CIV 2.176162 0.034715 NA
Ici, on voit que la Côte d’Ivoire est le seul pays ayant une influence sur le modèle. Mais la p-value est tout juste inférieure à 0.05 donc ce n’est pas très grave.